AD-A275  526 


r 


NASA  Contractor  Report  191583 
ICASE  Report  No.  93-95 


ICASE 


ON  THE  DAUBECHIES-BASED  WAVELET 
DIFFERENTIATION  MATRIX 


Leiand  Jamescn 


NASA  Contract  No.  NASl-19480 
December  1993 

Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  Virginia  23681-0001 

Opmted  by  the  Universities  Space  Research  Association 


National  Aeronautics  and 
Space  Administration 

Langley  Research  Center 

Hampton,  Virginia  23681  -0001 

94  2  l4  024 


94-04983 


ON  THE  DAUBECHIES-BASED  WAVELET 
DIFFERENTIATION  MATRIX 


Leland  Jameson} 

Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  VA  23681 

ABSTRACT 

The  differentiation  matrix  for  a  Daubechies-based  wavelet  basis  will  be  constructed 
and  ‘superconvergence’  will  be  proven.  That  is,  it  will  be  proven  that  under  the 
assumption  of  periodic  boundary  conditions  that  the  differentiation  matrix  is  accurate 
of  order  2M,  even  though  the  approximation  subspace  ran  represent  exactly  only 
polynomials  up  to  degree  M  —  1,  where  M  is  the  number  of  vanishing  moments  of  the 
associated  wavelet.  It  will  be  illustrated  that  Daubechies-based  wavelet  methods  are 
equivalent  to  finite  difference  methods  with  grid  refinement  in  regions  of  the  domain 
where  small-scale  structure  is  present. 


^This  research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under 
NASA  Contract  No.  NASl-19480  while  the  author  was  in  residence  at  the  Institute  for  Computer 
Applications  in  Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA 
23681.  Research  was  also  supported  by  AFOSR  grant  93-1-0090,  by  DARPA  grant  N00014-91-4016, 
and  by  NSF  grant  DMS-9211820. 


1 


1  Introduction 


The  term  differentiation  matrix  was  coined  by  E.  Tadmor  in  his  review  on  spectral 
methods  [1].  The  term  denotes  the  transformation  between  grid  point  values  of  a 
function  and  its  approximate  first  derivative.  This  matrix  is  a  product  of  three  ma¬ 
trices. 

The  first  matrix  C  is  constructed  as  follows:  assume  that  the  point  values  of  a 
function  /(x)  (where  a  <  x  <  b)  are  given  at  N  points  x,  for  0  <  j  <  N  —  1. 
Thus  a  vector  of  numbers  f(xj )  is  given.  FVom  this  vector  one  can  reconstruct  an 
approximation  to  the  fimction  /(x)  for  every  point  x  in  the  interval.  This  approxima¬ 
tion  (denoted  by  Psf)  itself  belongs  to  a  finite  dimensional  space  -  in  pseudospectral 
methods  it  is  the  global  interpolation  polynomial  that  collocates  f{xj)  and  in  finite 
differences  or  finite  elements  it  is  a  piecewise  polynomial.  This  transformation  be¬ 
tween  f{xj)  and  Pjv/j  defines  the  matrix  C.  Of  course  this  matrix  depends  on  the 
special  basis  chosen  to  represent  Psf^  A  good  example  is  the  Fourier  interpolation 
procedure  in  which  the  basis  is  the  set  of  complex  exponentials. 

The  second  matrix  D  results  from  differentiating  Psf,  and  projecting  it  b2w:k  to 
the  finite  dimensional  space.  Thus  D  is  defined  by  the  linear  transformation  between 
PnJ  and  Pn^Pn/- 

The  last  matrix  is  the  inverse  of  the  first  matrix.  Basically,  since  we  are  given  the 
approximation  we  can  read  it  at  the  grid  points  Xj.  Thus  the  differentiation 

matrix  T>  can  be  represented  as  I?  =  C~^DC. 

In  this  paper  the  wavelet  differentiation  matrix  will  be  examined.  As  with  other 
basis  sets,  as  outlined  above,  it  is  a  product  of  three  matrices.  Under  the  assumption 
of  periodicity  of  /(*),  however,  the  matrices  C  amd  D  commute  allowing  D  to  operate 
directly  on  the  vector  of  numbers  f{xj).  That  is,  the  differentiation  matrix  T>  is  simply 
D:  T>  =  D.  Furthermore,  the  matrix  D  differentiates  samples  of  polynomials  exactly, 
i.e.,  the  action  of  D  is  equivalent  to  a  finite  difference  operator  with  order  of  accuracy 
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d^Moidiiig  on  the  order  of  the  wevelet  choeen. 

More  precisely,  the  following  outlines  the  proof  of  this  assertion: 

•  Given  a  periodic  function  f{x),  let  C  be  the  mapping  from  evenly-spaced  sam¬ 
ples  of  f{x)  to  the  approximate  scaling  function  coefficients  on  the  finest  scale: 
C  :  f  (Tq.  Due  to  the  periodicity  of  /(x),  C  is  circulant  in  form. 

•  Let  D  be  the  mapping  from  the  exact  scaling  function  coefficients  of  /(x)  to 
scaling  function  coefficients  of  /'(x):  D  :  io  -♦  Sq.  Once  agun,  due  to  the 
periodicity  of  /(x),  D  is  circulant  in  form. 

•  The  matrix  operator  D  can  differentiate  exactly  evenly-spaced  samples  of  poly¬ 
nomials,  i.e.,  D  has  the  effect  of  a  finite-difference  operator.  The  order  of  exact 
differentiation  depends  on  the  order  of  the  wavelet  used. 

•  All  circulant  matrices  with  the  same  dimensions  conunute,  therefore  the  oper¬ 
ator  D  can  be  applied  directly  to  /: 

f  =  C-^DCf, 

or  simply, 

/'  =  Df, 

and  this  will  complete  the  proof. 

This  paper  contains  five  sections; 

§1)  This  introduction. 

§2)  Standard  definitions  of  wavelets  and  scaling  functions  are  given. 

§3)  The  general  approximation  properties  of  wavelets  will  be  discussed  along  with 
the  quadrature  formula  needed  to  approximate  the  scaling  function  coefficients  of  a 
function. 


|4)  This  is  the  most  important  section  of  this  paper.  It  will  be  proved  that  the 
action  of  D  is  equivalent  to  a  finite  difference  operator. 

§5)  The  results  of  sections  (3)  and  (4)  are  combined  for  the  desired  conclusion. 

In  addition,  the  foUowing  two  related  topics  are  explored  in  the  first  two  appen¬ 
dices: 

Appendix  A)  For  wavelets  supported  on  (0,3M)  it  will  be  shown  that  /  <l>{x)x"*dx  = 
(/  <^{x)xdx)”^. 

Appendix  B)  The  moments  of  the  scaling  function  ^(z)  will  be  calculated. 
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2  Wavelet  Definitions  and  Relations 

The  term  wavelet  is  to  describe  a  spatially  localized  function.  ‘Localized’  means 
that  the  wavelet  has  compact  support  or  that  the  wavelet  almost  has  compact  sup¬ 
port  in  the  sense  that  outside  of  some  interval  the  amplitude  of  the  wavelet  decays 
exponentially.  We  will  consider  only  wavelets  that  have  compact  support  and  that 
are  of  the  type  defined  by  Daubechies  [2]  which  are  supported  on  [0,2Af  —  1],  where 
M  is  the  number  of  vanishing  moments  defined  later  in  this  section. 

To  define  Daubechies  wavelets,  consider  the  two  functions  <^{x)  and  ij>{x)  which 
are  solutions  to  the  following  equations: 

^(i)  =  V2  £  ki,4(2x  -  k),  (1) 

J:=0 

v>(x)  =  \/2  9kH2x  -  k),  (2) 

k=0 

where  ^(z)  is  normalized, 

f  <l>{x)dx  =  1.  (3) 

00 

Let, 

^j(i)  =  -  k),  (4) 

and 

0i(x)  =  2"i^(2"^z  -  *),  (5) 

where  j,  k  ^  Z,  denote  the  dilations  and  translations  of  the  scaling  function  and  the 
wavelet. 

The  coefficients  H  —  and  G  —  are  related  by  gk  =  (— l)*hj[,_*  for 

ib  =  0, ...,  L  —  l.  Furthermore,  /f'and  G  are  chosen  so  that  dilations  and  translations 
of  the  wavelet,  ^(x),  form  an  orthonormal  basis  of  L^{R)  and  so  that  V’(^)  has  M 
vanishing  moments.  In  other  words,  0i(^)  satisfy 

(6) 

^—00 
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whwe  Su  is  the  Kronecker  delta  function.  Also,  rl>{x)  =  0o(^)  satisfies 

f  0(x)x"*<iz  =  0,  (7) 

y.oo 

for  m  =  0,  —  1.  Under  the  conditions  of  the  previous  two  equations,  for  any 

function  /(x)  €  L^(R)  there  exists  a  set  such  that 

/(*)  =  (8) 

j^ZKez 

where 

^jk  =  I  fix)ili{x)dx.  (9) 

•/— 00 

The  number  of  vanishing  moments  of  the  wavelet  tj){x)  defines  the  accuracy  of 
approximation.  The  two  sets  of  coefficients  H  and  G  are  known  in  signal  processing 
literature  as  quadrature  mirror  filters  [3].  For  Daubechies  wavelets  the  number  of 
coefficients  in  H  and  G,  or  the  length  of  the  filters  H  and  G,  denoted  by  L,  is  related 
to  the  number  of  vanishing  moments  M  by  2M  =  L.  For  example,  the  famous  Haar 
wavelet  is  found  by  defining  H  aa  ho  —  hi  =  1.  For  this  filter,  H,  the  solution  to 
the  dilation  equation  (1),  ^(x),  is  the  box  function:  ^(x)  =  1  for  x  €  [0, 1]  and 
^(x)  =  0  otherwise.  The  Haar  function  is  very  useful  as  a  learning  tool,  but  it 
is  not  very  useful  as  a  basis  fimction  on  which  to  expand  another  function  for  the 
important  reason  that  it  is  not  differentiable.  The  coefficients,  needed  to  define 
compactly  supported  wavelets  with  a  higher  degree  of  regularity  can  be  found  in  [2]. 
As  is  expected,  the  regularity  increases  with  the  support  of  the  wavelet.  The  usual 
notation  to  denote  a  Daubechies  wavelet  defined  by  coefficients  H  of  length  L  is  Di. 

It  is  usual  to  let  the  spaces  spanned  by  ^(x)  and  V^(x)  over  the  parameter  k, 
with  j  fixed,  to  be  denoted  by  Vj  and  Wj  respectively: 


The  speoet  Vj  end  Wj  ere  rdeted  by  [2] 


...  C  V,  C  Vo  C  V_1  C  ...»  (12) 

end  thet 


=  (13) 

The  previously  stated  condition  thet  the  wavelets  form  an  orthonormal  basis  of  L^{R) 
can  now  be  written  as, 

m  =  ©»Vi.  (14) 

i€Z 

Two  final  properties  of  the  spaces  Vj  are  that 


n  = w*  (15) 

jez 


and 


U  v^  =  L\R).  (16) 

i^z 

Properties  of  the  Semi-Discrete  Fourier  Transform  (SDFT)  of  the  filter  H  will  also 
be  needed.  The  following  definition  is  not  exactly  the  SDFT  but  a  constant  times  the 
SDFT: 

ibsL-l 

A({)  =  2-'/»  £  (IT) 

This  DFT  satisfies  the  following  equation,  see  [4]: 

\H{0\^  +  \H{(  +  ir)\^  =  l.  (18) 

S<dutions  of  equation  (18)  have  the  following  properties,  see  [2]: 

A(0  =  (5(1  +  «*))*'<?(«**),  (19) 

where  M  is  the  number  of  vanishing  moments  of  the  wavelet  and  Q  is  a  trigonometric 
polynmnial  such  that, 

W(«‘<P  =  i>(«n»{£/2))  +  «D“'((/2)Jt(l  co.(),  (20) 
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where 


and  A  is  an  odd  polynomial  such  that, 


0<P(y)  +  y*'i2(l/2-y) 


for  0  <  y  <  1,  and 


if  Af  >  2  or 


sup  (P(y)  +  y^R{ll2  —  y))  < 

0<y<l 


1-2H  l+2|x|’ 


for  |x|  <  1/2,  if  Af  =  1.  The  important  point  here  is  that  H{()  has  a  zero  of  order 
Af  at  ^  =  ir. 

Of  course,  infinite  sums  and  unions  are  meaningless  when  one  begins  to  implement 
a  wavelet  expansion  on  a  computer.  In  some  way  one  must  limit  range  of  the  scale 
parameter  j  and  the  location  parameter  k.  Consider  first  the  scale  parameter  j.  As 
stated  above,  the  wavelet  expansion  is  complete:  L^{R)  =■  Therefore,  any 

/(x)  €  L^{R)  can  be  written  as, 

/(®)  =  E  E 

i^Zk^Z 

where  due  to  orthonormality  of  the  wavelets  /(x)V’fc(x).  In  this  expan¬ 

sion,  functions  with  arbitrarily  small-scale  structures  can  be  represented.  In  practice, 
however,  there  is  a  limit  to  how  small  the  smallest  structure  can  be.  This  would 
depend,  for  example,  on  how  fine  the  grid  is  in  a  numerical  computation  scenario  or 
perhaps  what  the  sampling  frequency  is  in  a  signal  processing  scenario.  Therefore, 
on  a  computer  an  expansion  wotild  take  place  in  a  sp^tce  such  as 


Vq  =  Wx®W2®...®Wj®Vj, 
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and  would  ^pear  as, 


Pvofi^)  =  +  E  S  4v»fc(®),  (26) 

kez  j=i  kez 

where  again  due  to  orthonormality  of  the  basis  functions  /(x)V^(x),  and 

H  —  f^oo  ^  expansion,  scale  j  =  0  is  arbitrarily  chosen  as  the  finest 

scale  that  is  needed,  and  scale  J  would  be  the  scale  at  which  a  kind  of  local  average, 
provides  sufficient  large  scale  information.  In  language  that  is  likely  to  appeal 
to  the  electrical  engineer  it  can  be  said  that  (ff'lix)  represents  the  direct  current  portion 
of  a  signal  and  that  ^{{x)  represents  the  alternating  current  portion  of  a  signal  at, 
very  roughly,  frequency  j.  As  stated  above,  one  must  also  limit  the  range  of  the 
location  parameter  k.  In  this  paper  this  is  done  by  assuming  that  /(x)  is  a  periodic 
function.  The  periodicity  of  /(x)  induces  periodicity  on  all  wavelet  coefficients,  sj 
and 

This  completes  the  definition  of  wavelets.  The  next  section  will  discuss  function 
approximation  in  a  wavelet  basis. 


8 


3  Approximating  in  a  Wavelet  Basis 

Scaling  functions  and  wavelets  were  defined  in  the  previous  section.  The  goal  of  this 
section  is  to  find  the  coefficients  in  a  wavelet  expansion.  More  precisely,  the  scal¬ 
ing  function  coefiBcients  at  the  finest  scale,  sq,  will  be  approximated.  The  key  to 
this  approximation  is  the  matrix  C  which  maps  evenly-spaced  samples  of  a  periodic 
function  to  the  approximate  scaling  function  coefficients.  This  matrix  C  has  the  de¬ 
sirable  property  of  being  circulant  in  form  with  the  ramification  that  C  will  commute 
with  auiy  other  circulant  matrix,  particularly  the  derivative  matrix  72°,  the  subject  of 
section  (4).  An  example  of  C  is  given  at  the  end  of  section  (3.2). 

The  scenario  for  this  section  is  as  follows:  let  the  finest  scale  be  scale  j  =  0,  i.e., 
at  this  scale  all  relevant  small  scale  structures  in  the  function  have  been  captured  and 
represented.  One  seeks  an  expansion  of  a  function  /(x)  in  terms  of  ^2  in  the  space 
Vo.  With  the  projection  /Vo  defined  as  Pvo  :  Ki  such  an  expansion  has  the 

following  form: 

=  (27) 

k€Z 

where  due  to  the  orthonormality  of  the  basis  functions,  ^°(x)^°(x)dx,  the 

coefficients  s®  are  given  by 

4  =  [  fix)<f>ki^)dx.  (28) 

V— OO 


Once  the  s®  have  been  found  one  would  usually  then  find  the  scaling  function  and 
wavelet  function  coefficients  at  more  coarse  scales.  This  can  be  done  by  using  equation 
(29)  to  get  for  j  =  1,...,  J  and  by  using  equation  (30)  to  get  for  j  =  1,...,  J. 
These  equations  are  derived  respectively  from  equations  (1)  and  (2),  see  [4],  [5], 

3M 

4  =  E  (29) 

n=l 


and 


2Af 

=  ^29n4+2k-2' 

n=l 


(30) 
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In  this  section,  however,  the  decomposition  onto  more  coarse  scales  will  not  be  cal¬ 
culated.  The  important  step  for  this  section  is  the  approximation  of  the  integral 
^2  =  /-^  ^2  the  approximation  to  s®.  The  quadrature  for¬ 

mula  for  this  integral  encompasses  the  approximation  properties  of  scaling  functions, 
and  hence  wavelets. 

This  section  contains  3  subsections: 

§3.1  )  The  approximation  properties  of  scaling  functions  will  be  discussed  and  the 
quadrature  formula  to  estimate  the  integral  f{x)<f>^{x)dx  will  be  derived. 

§3.2  )  An  example  using  the  results  from  section  (3.1)  is  given  for  the  Daubechies 
wavelet  De- 

§3.3  )  The  example  from  section  (3.2)  leads  to  a  circulant  matrix  for  the  matrix 
C,  Circulant  matrices  will  be  defined  and  the  ramifications  of  circularity  will  be 
discussed. 

3.1  Quadrature  Formula  for  Scaling  Function 

In  this  subsection  the  coefiBcients  s*  will  be  approximated.  Before  stating  the  appro¬ 
priate  quadrature  formula,  however,  the  order  of  accuracy  of  a  wavelet  approximation 
is  discussed. 

3.1.1  Approximation  Properties  of  Scaling  iXinctions 

This  subsection  comes  essentially  from  [6].  The  approximation  properties  of  scaling 
functions  are  determined  by  the  Discrete  Fourier  Transform  of  the  filter  H.  That  is, 
if 

=  (3>) 

^  *=o 

has  a  zero  of  order  Af  at  (  =  x  then  there  are  a  number  of  consequences: 

1.  The  polynomials  1,  x,  ...,  x^"'  are  linear  combinations  of  the  translates  of  the 
scaling  function 
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2.  Smooth  functions  can  be  approximated  with  error  0{h^),  where  h  denotes  the 
grid  size,  i.e.,  there  exist  a  set  where  j  is  fixed,  and  there  exists  a  constant 
C  such  that 

l/(*)  -  (32) 

k 

where  the  norm  i  ■  { is  the  Lj  norm. 

3.  The  associated  wavelet  has  M  vanishing  moments, 

J  x"^^{x)dx  =  0 

for  m  =  0,  —  1. 

Other  ramifications  can  be  found  in  [6].  These  approximation  properties  determine 
the  accuracy  of  the  quadrature  formula  used  to  approximate  the  scaling  function 
coefficients  io  which  is  derived  in  the  following  section. 

3.1.2  Derivation  of  Quadrature  Formula 

It  is  important  to  note  that  all  wavelets  in  this  paper  are  of  the  usual  Daubechies 
type,  i.e.,  the  support  of  a  usual  Daubechies  wavelet  D2M  is  [0,2Af  —  1]  where  M  is 
the  number  of  vanishing  moments  of  the  wavelet.  For  this  subsection  this  support  size 
is  partictilarly  important  to  keep  in  mind  because  there  does  exist  an  orthonormal 
family  of  wavelets  which  are  supported  on  [0, 3Af  —  1]  and  which  have  a  very  simple 
quadrature  formula  based  on  the  vanishing  moments  of  the  wavelet  (see  appendix  A) 
but  this  is  not  the  wavelet  being  used  in  this  paper. 

Given  the  approximation  properties  of  the  scaling  function  from  the  previous 
subsection,  one  can  now  seek  a  quadrature  formula  which  is  exact  when  /(x)  is  a 
polynomial  up  to  order  M  —  1:  f{x)=  p(x)  €  Pju-i-  That  is,  there  exist  a  set  of 
coefficients  {c/}j^o  ^  such  that 

/  P{x)<i>ldx  =  ^  cip{l  +  k),  (33) 

/=o 
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for  p(x)  €  Pm-1.  If  the  integral  is  shifted  the  above  equation  becomes, 

/  P{y  +  k)<f>l{y)dy  =  ^  cjp(/  +  A:).  (34) 

•'-«»  {sO 

More  simply,  the  coeifficients  can  be  found  [9]  by  solving  the  following  linear 

system: 

«!)&  =  ^  rc,  (35) 

issO 

for  m  =  -  1,  and  the  coefficients  {c/}j^^  provide  the  desired  quadrature 

formula.  That  is,  the  coefficients,  cq,  which  approximate  sq  are  found  from, 

<’1  =  E  <=>/(' + *)•  W 

1=0 

When  placed  in  matrix  form  the  coefficients  {cj}j^*  yield  the  circulant  matrix  C. 
A  more  thorotigh  discussion  of  circulant  matrices  will  be  given  in  §(3.3)  after  the 
example  of  the  next  subsection  has  been  completed. 

Note  that  since  the  above  derived  quadrature  formula  is  exact  for  p(x)  6  Pm-i 
the  coefficients  <t®  approximate  the  coefficients  s®  with  error  of  order  M .  Also,  note 
that  the  derivation  of  the  quadrature  coefficients  depends  only  on  the  moments  of  the 
scaling  function,  x"*^(x)dx.  In  the  next  subsection,  the  moments  of  the  scaling 
function  will  first  be  calculated  and  then  the  coefficients  will  be  found  for 

the  Z?e  wavelet.  The  wavelet  De  is  chosen  for  no  other  reason  tham  that  D2  and  D4 
receive  considerable  attention  from  other  sources  and  that  Dg  is  slightly  less  trivial 
than  D2  and  D4  while  remaining  manageable. 

3.2  Example  with  the  Daubechies  Wavelet 

Recall  from  the  previous  subsection  that  the  inunediate  goal  is  to  approximate  the 
scaling  function  coefficients  of  4  function  at  scale  j  =  0.  Specifically,  in  this  section 
the  objective  is  to  derive  the  matrix  form  of  the  mapping  from  evenly-spaced  samples 
of  a  periodic  fimction  /(x)  to  the  scaling  function  coefficients  on  the  finest  scale  s®. 
The  example  will  be  for  the  Daubechies  wavelet  Dg.  Comparable  results  for  the 
wavelets  D4  and  Ds  are  presented  in  appendix  B. 
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Rjecall  from  the  previous  subs^tion  that  in  order  to  calculate  the  coefiSdents 
the  moments  of  the  scaling  function  ^(x)  must  first  be  known.  Let  Mi  be 
the  moment  of  the  scaling  function  ^x), 

M,  =  J  ^(x)x'dx,  (37) 

and  let  m  be  the  moment  of  the  filter 

(38) 

k 

The  zero-th  moment,  Afo,  of  ^(x)  is  1  by  the  normalization  of  ^(x): 


Mq  =  j  ^(x)dx  =  1. 


The  zero-th  moment  of  the  coefficients  hk  is  fotmd  by  integrating  the  dilation  equation 


which  defines  ^(x): 


J  <f>{x)dx  =  j  ^(2x  —  fc)dx, 


and  let  y  =  2x  —  A;  to  get, 


which  implies, 


=  jE**/  miy. 


(<o  =  E*“  =  2-  («) 

k 

That  is,  the  zero-th  moments  Mq  and  fiQ  are  the  same  for  all  Daubechies  wavelets. 
Higher  moments  for  /ij  can  be  found  by  straight-forward  calculation  using  the  coeffi¬ 
cients  provided  by  Daubechies  [2].  The  higher  moments,  Mi  for  1  >  0,  for  the  scaling 
function  can  be  found  from  the  following  equation  which  is  derived  in  appendix  B: 


=  (j)”-"'  E  ( 7 ) 


For  the  current  example  only  the  moments  A/o,  M\,  and  M^  are  needed;  Mq  =  1, 
Ml  =  5^1,  and  Afj  =  g((/ii)^  +  ^2)-  Given  these  three  moments  the  coefficients 
{c/}j^o^  can  be  fotmd  from 


^  r*c,  =  /x”*^(x)dx 

f_A  J 
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i  Mi  fij  Cj 

0  1  2  .1080 

1  .8174  1.6348  .9667 

2  .6681  1.3363  -.0746 

Table  1:  Scaling  function  and  filtei  moments  for  Daubediies  6  wavelets. 

for  m  =  0, 1, ...,  M  —  1.  Specifically,  for  the  Dg  wavelet  the  linear  system  in  matrix 
form  is, 

(;  I fs). 

\0  1  4)  \  c2)  {Mil 


which  has  the  solution  cq  =  .1080,  ci  =  .9667,  and  cj  =  —.0746.  In  tabular  form,  the 
complete  results  for  Da  are. 

Recall  that  the  quadrature  formula  used  to  approximate  the  scaling  function  has 
the  form, 

<’t  =  E'  Cif(l  +  *)•  («) 

1=0 


If  the  function  f(x)  is  periodic  then  in  matrix  notation  the  above  operation  is  =  (7/ 
where  C  for  Da  and  on  a  grid  of  6  points  is, 


.108 

.967 

-.075 

0 

0 

0 

0 

.108 

.967 

-.075 

0 

0 

0 

0 

.108 

.967 

-.075 

0 

0 

0 

0 

.108 

.967 

-.075 

-.075 

0 

0 

0 

.108 

.967 

.967 

-.075 

0 

0 

0 

.108 

The  important  point  here  is  that  the  above  matrix  is  circulant.  The  ramifications 
of  circularity  are  very  important  for  the  thesis  of  the  paper.  The  definition  of  circulant 
matrices  and  the  properties  that  they  are  imbued  with  is  the  subject  of  the  next 
subsection. 

3.3  Circulant  Matrices 

Strang  [7]  defines  a  circulant  matrix  as  a  constant-diagonal  matrix  which  is  ‘periodic, 
since  the  lower  diagonals  fold  around  to  appear  again  as  the  upper  diagonals.’  A 
thorough  discussion  of  circulant  matrices  is  given  by  Davis  [8].  Circulant  matrices 
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have  the  wcmderful  property  that  they  can  all  be  diagonalized  by  the  same  matrix,  the 
Fourier  matrix:  An  NxN  Fourier  matrix  has  as  its  ij —tk  element  the  entry 
where  vo^  —  1.  The  most  important  ramification  for  this  paper  is  that  matrices  which 
can  be  diagonalized  by  the  same  matrix  conunute.  That  is,  the  matrix  C  from  the 
previous  subsection  will  commute  with  the  matrix  /2°  which  will  be  derived  in  section 
(1.4). 

In  general,  circulant  matrices  arise  whenever  one  is  performing  the  matrix  version 
of  periodic  discrete  convolution.  In  munerical  analysis  periodic  discrete  convolution 
arises  whenever  one  differentiates  the  evenly-spaced  samples  of  a  periodic  equation 
which  has  constant  coefficients.  Let  us  be  a  bit  more  precise  and  illustrate  how  the 
operation  of  periodic  discrete  convolution  yields  a  circulant  matrix  by  stating  the 
following  theorem: 

Theorem:  A  finite-length  filter  of  length  M  applied  to  N  evenly-spaced  samples 
of  periodic  fimction,  where  N  >  M,  will  in  matrix  form  yield  a  circulant  matrix. 

Proof:  First  of  all,  let  the  notation  remain  as  above:  co,  ci, ...,  cm-i  will  represent 
the  finite-length  filter  and  /o,  /i, ...,  fpt-i  will  represent  the  evenly-spaced  samples  of 
one  period  of  the  periodic  function  /(x).  Of  course,  the  samples  of  /(x)  are  periodic 
with  period  N.  The  application  of  the  filter  c  on  the  samples  of  /(x)  is  the  convolution: 

M-l 

<^k=  ^ifk-h 
/=0 

A 

where  /<  is  the  renaming  of  the  elements  of  /<  so  that  the  previous  convolution  is 

A 

the  same  as  the  following  expression.  Furthermore,  keep  in  mind  that  fi  and  fi  are 
periodic  with  period  N. 

O’*  =  53 
1=0 

Using  the  modulus  function  to  keep  the  indices  of  /,  within  one  period,  i.e.,  keep 
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0  <  t  <  —  1,  the  above  equation  can  be  written  as, 


M-t 

—  ^  Clfmo^t+kjny 
ImO 

Now,  shift  the  indices  by  letting  j  =i  I +  k  to  get, 

M-(M-l) 

j** 

If  the  length-M  filter  c  is  now  ’padded*  at  the  end  with  zeros  so  that  it  is  now  a 
length-N  filter  then  the  above  equation  can  be  rewritten  as, 

s-\ 

i=o 

This  is  exactly  a  matrix  multiply  S  —  Cf  where  the  ij  —  th  element  of  the  matrix 
C  is  wd  this  is  the  definition  of  a  circulant  matrix.  This  completes  the 

proof./ / 

Note  that  the  difference  between  a  circulant  matrix  and  a  Toeplitz  matrix  is  the 
wrapping  arotmd  effect  of  the  diagonals  introduced  by  the  use  of  the  modulus  function 
for  the  circulant  matrix.  That  is,  a  circulant  matrix  is  a  special  case  of  a  Toeplitz 
matrix  where  the  constant  diagonals  are  periodic. 

Before  leaving  the  discussion  on  circulant  matrices  let  one  more  interpretation  be 
noted:  to  say  that  circulant  matrices  conunute  is  to  simply  restate  the  important  re¬ 
sult  from  signal  analysis  that  convolutions  conunute.  That  is,  if  one  has  two  sequences 
c  and  r,  which  in  the  current  scenario  are  periodic,  then  the  order  of  convolution  does 
not  matter.  This  is  easily  proved  with  the  Fourier  transform: 

o'*" r  =  cf  =  fc  =  r“*c,  (48) 

where  f  denotes  the  Fourier  transform  of  r  and  denotes  periodic  convolution.  For 
this  paper,  c,  of  course,  would  be  the  quadrature  operator  and  r  would  be  the  scaling 
frmction  derivative  operator  which  is  the  subject  of  section  (4).  In  matrix  notation 
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equation  (48)  is  nothing  more  than, 


C  •  /i“  =  •  C  (49) 

In  this  section  a  quadrature  formula  has  been  found  to  approximate  the  scaling 
function  coefficients  of  a  given  function,  f{x).  In  matrix  form  this  quadrature  for¬ 
mula  leads  to  a  circulant  matrix  assuming  f{x)  is  periodic.  In  the  next  section  the 
wavelet  derivative  operator  will  be  derived,  and  it  will  be  shown  that,  once  again, 
the  assumption  that  f{x)  is  periodic  leads  to  an  opo-ator  which  in  matrix  form  is 
circulant. 
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4  Derivative  based  on  Wavelets 

In  the  previous  section  the  mapping  from  evenly-spaced  sunples  of  a  periodic  function, 
/(x),  to  the  i^proximate  scaling  function  coefficients  on  the  finest  scale,  (r°,  was  given. 
Recall  that  denotes  the  apprmcimation  to  the  exact  scaling  function  coefficient 
at  scale  j  and  position  k.  The  mi^ping  is  nothing  more  than  a  quadrature  formula 
which  is  exact  when  f{x)  is  equal  to  a  polynomial  up  to  order  Af  —  1,  where  M  is 
the  number  of  vanishing  moments  of  the  wavelet.  The  question  now  is  how  does  one 
represent  the  derivative  of  /(x)  in  the  wavelet  basis  ^ven  that  the  wavelet  expansion 
of  /(x)  is  already  given. 

The  answer  is  given  in  the  following  subsections: 

§4.1  )  A  function  /(x)  will  be  expanded  in  a  wavelet  basis  and  the  expansion  will 
be  differentiated. 

§4.2  )  The  results  from  Beylkin  [9]  on  derivative  projections  will  be  ^ven. 

§4.3  )  First,  it  will  be  noted  that  one  can  differentiate  a  wavelet  expansion  at  any 
level  of  a  wavelet  decomposition  and  achieve  the  same  derivative.  Second,  explicit 
wavelet  decomposition  will  be  performed  accompanied  by  the  appropriate  differenti¬ 
ation  matrix. 

§4.4  )  The  similarity  between  the  wavelet-based  derivative  coefficients  and  finite 
difference  derivative  coefficients  will  be  noted,  and  it  will  be  shown  that  when  one 
differentiates  the  wavelet  expansion  of  a  periodic  function  that  the  effect  on  the 
original  function  samples  is  equal  to  finite  difference  differentiation. 

4.1  Expansion  in  a  Wavelet  Basis 

The  goal  now  is  to  find  the  wavelet  and  scaling  function  expansion  of  a  periodic 
function  /(x).  Given  /(x)  €  L^{R)  one  first  projects  onto  the  arbitrarily  chosen 
finest  scale  j  =  0  of  the  scaling  function  <f>°{x)  which  generates  the  space  Vo,  i.e.,  let 
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(50) 


be  the  i»rojectkm  frmn  the  space  L^{R)  to  the  space  Vo,  /Vo  :  -*  Vo- 

A4/(*)  =  E  *!♦!(*). 

kaa/Q 

where  due  to  the  orthonormality  of  ^ 

«*  =  /  (*)<*»•  (51) 

Note  that  in  the  introduction  the  projection  denoted  by  Pn  would  be  the  same  as 
Pv^  using  notation  that  is  more  amenable  to  wavelets.  The  derivative  of  Pv^f{x)  is, 

^fl4/(*)  =  E  <2^S(*)-  (“) 

Of  course,  iPvJ{x)  is  not  in  Vq  and  must  be  projected  onto  Vq.  First  define  the 
inner  product  <  f^g>  <m  L^(R)  by 

<f,9>-f  (53) 

•r— 00 

Now  the  projection  of  ^/Vo/(*)  onto  Vo  is, 

/V.^/V./(x)=£<^iV./,1>f>^?(*),  (54) 

IsO 

or, 

ft..^ft4/(*)=E  E*J<^2.*?>*f(»)-  (55) 

(sO  *s0 

_  # 

The  inner  product  <  ^2,  >  is  one  of  the  results  provided  in  [10]. 

In  the  previous  paragri^h  /(x)  was  expanded  in  a  scaling-function  expansion  at 
the  finest  scale  j  =  0.  Now  /(x)  will  be  expanded  in  terms  of  scaling  functions  and 
wavelets  at  scale  j  =  1 .  Recall  that  =  V^  0  IFi.  Now  one  must  project  from  L^{R) 
onto  V^  and  from  L^{R)  onto  W\.  Let  both  projections  be  denoted  simultaneously 
by  That  is,  /View,  :  L^R)  -*  Vi©W'i.  Let  /V,ewi/(*)  be  the  projection 

of  /(x)  on  V^  0H^i.  Therefore,  the  expansion  for  /V,ew,/(x)  is, 

N/2-X  N/2-1 

A'..m/(*)=  E  *!«(*)+  E  (56) 

kaO  k=0 
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iribm  <lne  to  th«  orthaaocmolity  oC  the  bast  fiuictkxu  and  ^k(*)  ^  oo^deota 

•I  and  ^  am  pna  by 

•k^  r  /(*¥i(*)d», 

«/— OP 

and 

iT— 00 

The  dcrivetive  of  if 


(57) 


(58) 


.  w/i-i  ^  w/a-i 

■^PviBWifix)^  2  44(*)+  S 

**  "  k«0 


(59) 


Once  a^n,  the  deriyative  of  Pvi9Wif{x)  does  not  belong  to  Vi0)Vi,  and  must, 
timrefom,  be  projected  back  onto  this  space.  The  projection  is, 

Pvi^Wi  ^/Viom  /(*)  =  (60) 

N/a-iN/a-i 

laO  farf) 

JV/a-iJV/a-i 

+  S  E  >W(a:) 

/am  AaiO 

jv/a-1  JV/a-1 

laO  SaO 

N/a-1  N/a-1 

+  E  E  >^?(x). 

(sO  kacO 

The  four  inner  products  <  >,  <  >»  <  >»  ®“*d  <  >  are 

the  key  to  finding  the  derivative  of  a  wavelet  expansion,  and  are  provided  in  [9].  An 
outline  of  the  derivation  of  these  inner  products  is  ^ven  in  the  next  section. 


4.2  Wavelet  Coefficients  of  the  Derivative 


An  arlntrary  wavelet  expansion  of  a  function  might  contain  wavelet  coeffidents  and 
scaling  coeffidents  at  many  scales.  In  [9]  the  projection  coefficients  that  map  from 
scaling  function  coeffidents  and  wavelet  function  coefficients  at  a  (pven  scale  to  the 
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derivative  scaling  function  coefficients  and  wavelet  function  coefficients  at  the  same 
scale  are  derived.  The  matrix  elements  of  these  projections  are  computed  from, 


2~^o,_{  =  aji  =  2“^^  f  xl}{2~^x  —  t)^(2~^x  —  l)dx, 

•/— OO 

(61) 

2-^bi.i  =  bi,  =  2-^-'  r  ^(2--'x  -  -  /)dx, 

J—oo 

(62) 

2-^Ci.i  =  4  =  2-^^  r  ^(2-^x  -  i)ti>{2-^x  -  l)dx, 

•/— OO 

(63) 

/oo  ^ 

=  4  =  2-^^  /  ^(2--'x  -  0^(2"-'x  -  l)dx. 

J—OO 

(64) 

Since  the  above  projections  2ire  always  at  a  fixed  scale,  j,  the  projection  coefficients 

are  simply. 

(65) 

(66) 

(67) 

“  /  oo  ~ 

(68) 

for  /  €  2. 

Furthermore,  using  the  dilation  equation  which  defines  ^(x),  ^(x)  = 

Y^k  hk<f>{2x  ■ 

-  A:),  and  the  equation  which  defines  rl^{x),  ^{x)  =  Yi,k9k<f>{2x 

—  k),  the 

first  three  of  the  above  four  equations  become. 

L~\  L-l 

<*«  =  5^  13  9k9ir2i+k-l 
kzsO  <=o 

(69) 

L-l L-l 

=  X/  53  9kbir2i+k-i 

(70) 

L-l  L-l 

Ct  =  13  53  f^k9l^2i-kk-l- 

fc=0  {=0 

(71) 

It  is  apparent  from  the  above  equations  that  the  coefficients  rj  contain  all  the  infor¬ 
mation  concerning  the  derivative.  The  coefficients  rj  can  be  found  [9]  from  solving 
the  following  system  of  linear  algebraic  equations: 


j  I</2 

n  =  2(r2/  +  X  ®2fc-i(r2j_2fc+i  -I-  r2/+2fc_i)),  (72) 

^  fc=i 
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*h1 


-  -1,  (73) 

< 

where 

a,  =  2  'Z  (74) 

tsO 

for  n  =:  1, ...,  X  —  1.  The  proof  of  the  above  pr<^>oeition  can  be  found  in  [9]. 

This  section  has  pven  a  brief  outline  of  the  derivation  of  the  wavelet  derivative 
projection  coefficients.  It  is  impcfftant  to  note  that  all  the  information  for  the  wavelet 
derivative  is  contained  in  the  coefficients  {ri},  and  this  point  will  be  explored  more 
in  next  section. 

4.3  Derivative  at  Scale  Zero  of  Scaling  F^mirtion  Only 

Wavelet  derivatives  can  be  calculated  at  any  level  of  a  wavdet  decomposition.  The 
result  will,  of  course,  always  be  the  same.  That  is,  recall  the  leUtioo  from  §(2), 
Vj  ss  Vj+i  .  As  stated  before,  it  is  the  conv«iti<ni  of  this  paper  to  let  Vo  represent 
the  finest  scale.  Using  the  above  relation,  (me  could  (fecompoee  any  number  of 
times.  One  decomposition  yields  V(i  =  Wi  0  Vi,  and  a  second  decomp(Miti<m  yields 
=  Wi  0  Wj  0  Va.  One  could  calculate  the  wavelet  derivative  in  any  one  of  these 
spaces.  Once  again,  the  goal  of  this  paper  is  to  understand  the  essence  of  a  wavelet 
derivative,  and  since  the  derivative  is  the  same  regardless  of  the  decomposition  of  the 
space,  one  should  ch(X)se  the  simplest  approach  and  calculate  the  derivative  at  scale 
j  =  0  using  only  the  seeding  function  coefficients. 

This  subsection  contains  four  parts: 

1.  New  notation  will  be  introduced. 

2.  Wavelet  decompositions  and  differentiation  matrices  will  be  given  for  the  space 
V{|  as  well  as  conunents  on  data  compression  in  this  space. 

3.  Wavelet  decompositions  and  differentiation  matrices  will  be  given  for  the  space 
Wi  0  Vi  as  well  as  comments  on  data  compression  in  this  space. 
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4.  Wavelet  decompoeitions  and  differentiation  matrices  will  be  given  for  the  space 
Wi  0  W2  0  W3  0  14  as  well  as  comments  on  data  compression. 


4.3.1  New  Notation 


To  simplify  the  presentation,  matrix  notation  will  be  used  whenever  possible.  Begin 
by  defining  the  matrix  version  of  equations  (29)  and  (30).  RecaU  that  these  equations 
are 


n=iM 

nasi 


and 


~  ^  9n^+2k-2' 

n=l 

Denote  the  decomposition  matrix  embodied  by  these  two  equations  by  where 

the  matrix  subscripts  denote  the  size  of  the  matrix,  and  the  superscripts  indicate 
that  P  is  decomposing  from  scaling  function  coefBcients  at  scale  j  to  scaling  function 
and  wavelet  function  coefficients  at  scale  j  +  1.  As  before,  let  Sj  contain  the  scaling 
function  coefficients  at  scale  j.  (Note:  When  vector  notation  is  used  the  scale  is 
given  as  a  subscript,  otherwise  the  location  k  is  the  subscript  and  the  scale  is  the 
superscript.)  P  therefore  maps  sj  onto  Sj+i  2Uid  dj+i". 


(75) 


Note  that  the  vectors  at  scale  j  +  1  are  half  as  long  as  the  vectors  as  scale  j.  To 
illustrate  further,  suppose  the  wavelet  being  used  is  the  four  coefficient  D4  wavelet, 
and  suppose  one  wants  to  project  from  8  scaling  function  coefficients  at  scale  j  to  4 
scaling  function  coefficients  at  scale  j  +  1  and  4  wavelet  coefficients  at  scale  j  +  1. 
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The  <ieo(»npoati<m  matrix  in  this  case  is, 

h]  h]  hs  ^4  0  0  0  0 

0  0  hi  h)  hs  ^4  0  0 

0  0  0  0  hi  h]  hs  ^4 

pjj+i  _  hs  h4  0  0  0  0  hi  hj  ,  . 

"  9i  93  94  0  0  0  0  • 

0  0  91  Sh  93  94  0  0 

0  0  0  0  gi  9i  9i  §4 

93  94  0  0  0  0  9i  ih  ^ 

Other  decomposition  matrices  of  different  sizes  will  have  the  same  structure  as  the 
above  matrix. 

For  a  bit  more  matrix  notation,  let  the  four  matrices  i4j^xN>  ^NxNi  ^NxNt 
I^NxN  contain  the  derivative  projection  coefficients  defined  in  §(4.2)  where,  again, 
the  subscripts  denote  the  size  of  the  matrix  and  the  superscript  denotes  the  scale  on 
which  the  derivative  projection  is  acting.  The  elements  of  the  four  matrices  are, 

A  Cij  =  Oj.j, 

B  ^  bfj  — 

C  c%j  —  c,_j. 


R  m  =  ri_„ 

and  the  mappings  performed  by  the  matrices  are. 


:  dj  dj, 

:  Sj  — ♦  dj, 

:  dj  -*  Sj, 

:  Sj  -*  ij  , 

where  Sj  and  d,-,  as  before,  define  the  scaling  and  wavelet  coefficients  at  scale  j,  and 
Sj  and  dj  denote  the  coefficients  of  the  expansion  of  the  derivative  of  a  function  which 
is  initially  defined  by  an  expansion  in  Sj  and  dj. 

This  concludes  the  new  notation. 
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4.3.2  Wavelet  Expansion  and  Deri^^tive  in  Vo 

As  stated  previously,  one  can  calculate  the  derivative  of  a  wavelet  expansion  at  any 
level  in  the  wavelet  decomposition.  This  subsection  will  explore  the  first  of  three  of  the 
options.  To  be  explicit,  suppose  that  a  periodic  function  /(x)  has  been  approximated 
on  a  grid  with  16  scaling  function  coefficients  to  get  ioj  a.nd  for  the  current  argument 
assume  that  the  coefficients  have  been  calculated  ex2w:tly,  i.e.,  the  notation  io  will  be 
used  instead  of  di.  Furthermore  due  to  the  periodicity  of  /(x)  the  coefficients  so  will 
also  be  periodic.  The  coefficients  of  the  expansion  of  ^/(af)  in  Vo  are  found  from  so 
by  an  application  of  the  previously  defined  matrix  B^exis' 


which  completely  defines  the  derivative  of  /(x)  in  the  scaling  function  basis  at  scale 
i  =  0. 

For  data  compression  purposes,  the  space  Vo  is  not  a  good  space  to  work  in.  That 
is,  the  coeflScients  Jq  represent  the  equivalent  of  a  local  averages.  In  a  wavelet  basis, 
it  is  often  true  that  the  coefficients  of  local  high-frequency  oscillations  are  small  and 
can  be  set  to  zero  without  altering  the  character  of  the  function  being  represented, 
but  the  coeflBcients  of  local  averages  usually  represent  essential  information. 
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4.S.S  Wavelet  Expansion  and  Derivative  in  W\  0  Vi 

Consider  now  a  decomposition  of  the  vector  of  scaling  function  coefficients  sq  onto 
the  scaling  ftmction  and  wavelet  coefficients  at  scale  j  =  1  by  an  application  of  the 
matrix  PiVxie: 


As  before,  we  have  16  basis  functions  in  our  space  which  is  now  Vi  0  Wi  rather  than 
Vo-  In  order  to  calculate  the  coefficients  of  the  derivative  expansion  in  Vi  0  Wi  the 
following  projections  are  calculated: 

=  ^XS  ■  +  ^8x8  ■  (79) 


and 

•  di  +  Bgxg  •  si,  (80) 

where  A,  B,  C,  and  R  were  all  defined  in  the  previous  subsection.  A  more  concise 
way  to  represent  the  derivative  projections  is  in  matrix  notation: 

£l  _  ilgx8  ^8X8  .  ^ 

di  .  ^8X8  -^8X8  .  .  . 

If  one  now  applies  the  matrix  (Pjgxie)^  {T  denotes  transpose  and  hence  inverse  for 
this  unitary  matrix)  to  the  derivative  coefficients  at  scale  j  =  I  then  one  gets  the 
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derivative  coefficients  at  scale  j  =  0: 


(82) 


and  one  gets  exactly  the  same  coefficients  as  before  when  the  matrix 
applied  to  sq.  To  emphasize,  the  derivative  calculated  at  scale  j  =  0  and  the  derivative 
calculated  at  j  =  1  yield  exactly  the  same  result.  The  importance  of  this  observation 
is  that  in  order  to  understand  the  essence  of  the  wavelet  derivative  one  need  only  be 
concerned  with  the  action  of  the  matrix  lUj^xN  vector  Sq. 

For  data  compression  the  space  0  Vi  is  a  fur  space  to  work  in.  The  coefficients 
Si  of  the  basis  fimctions  in  V\  represent  local  averatges  just  as  the  coefficients  of 
the  basis  functions  in  the  space  Vo  do.  However,  the  basis  functions  in  Vi  have 
broader  support  than  the  basis  functions  in  Vo  and  therefore  represent  averages  over 
a  larger  amount  of  data  (twice  as  much  data  to  be  exact).  Therefore,  once  again 
the  coefficients  Si  usually  carry  essential  information.  The  coefficients  di  of  the  basis 
fimctions  in  the  space  W\^  on  the  other  hand,  carry  information  concerning  local 
oscillations.  That  is,  if  the  function  being  represented,  /(x),  is  globally  smooth  then 
the  coefficients  d\  will  be  near  zero  and  can  be  set  exactly  to  zero  without  altering  the 
character  of  /(x).  In  the  solution  of  nonlinear  partial  differential  equations  where  a 
sharp  gradient,  or  shock,  can  develop,  the  coefficients  d\  away  from  the  shock  would 
be  close  to  zero  whereas  the  coefficients  near  the  shock  would  be  large.  Therefore, 
representing  a  function  in  IVi  0  Vi  is  more  versatile  than  simply  staying  in  the  space 
V^.  Versatility  continues  to  be  enhanced  as  one  decomposes  into  more  and  more 
wavelet  subspaces  as  in  the  next  and  final  scenario. 


4.3.4  Wavelet  Expansion  and  Derivative  in  Hi  0  Wj  0  W3  0  V3 

Up  to  now  our  basis  functions  have  all  been  at  the  S2mie  scale,  i.e.,  initially  our 
basis  functions  were  contained  in  Vq,  and  in  the  second  scenario  the  basis  functions 
were  contained  in  Vi  and  Wi.  In  this  subsection,  however,  the  basis  functions  will 
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be  contained  in  spaces  at  three  different  scales:  Wi,  and  V3.  The  full  set 

of  co<^cients  in  this  case  and  all  the  appropriate  decompositions  leading  to  these 
coefficients  are, 


1 
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■  Si  ■ 

[41 

S* 

si 

si 

si 

si 

4 

si 

S4 

.  si  . 

l4j 

si 

si 

[41 

[41 

si 

si 

4 

4 

s? 

s\ 

4 

4 

si 

si 

ptA 

[4  J 

*4)U 

4 

si 

A' 

[41 

[4] 

sio 

di 

4 

4 

^11 

4 

4 

4 

4 

A 

4 

4 

Si4 

4 

4 

4 

sis 

4 

4 

4 

.  sfe  . 

[41 

L4J 

L4J 

In  matrix  form  the  projection  onto  the  coefficients  of  the  derivative  of  the  expansion 
is,  where  the  matrix  will  be  labeled  M, 


^X2  ^3x2  p3|3  /nt2 

d3  j3  ^4x4'^4x4 

"3x3  -^3X3 

M  = 

Bl^4iP4^4f  Al^, 

0 

I 


) 


pl.2 

*^8x8 


^8x8 


■^8x8 


(84) 
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and  M  p«rf<»ina  the  fcdlowing  mapping: 


For  data  compression,  this  is  the  most  useful  set  of  subspaces.  The  space  now  is 
represented  as  Wi  ®  W2  B  W3  ®  V3.  For  the  same  reasons  as  before  the  coefficients 
of  basis  functions  in  the  subspace  V3  cannot  be  ignored.  It  is  likely,  however,  that 
the  function  f{x)  being  represented  is  smooth  in  most  of  the  domain  allowing  one 
to  disregard  the  majority  of  the  coefficients  of  the  basis  functions  in  the  subspace 
WxBW^B  W3.  In  fact,  it  is  more  likely  that  the  coefficients  for  the  basis  functions 
in  Wi  will  be  negligible  than  for  the  coefficients  for  the  basis  functions  in  W3.  This 
is  because  the  basis  functions  in  W3  have  larger  support  than  the  basis  functions  in 
Wi  and  Wi. 

In  summary,  an  attempt  has  been  made  to  illustrate  that  the  derivative  coefficients 
of  a  scaling  and  wavelet  expansion  can  be  calculated  at  any  scale.  The  proper  set  of 
wavelet  subspaces  depends  on  the  problem  at  hand.  The  goal  for  this  author  is  to 
understand  exactly  what  wavelets  are  and  what  they  are  doing,  therefore,  scale  j  =  0, 
i.e.,  the  space  Vo,  provides  the  clearest  scenario  in  which  to  work  without  sacrificing 
essential  properties  of  wavelets. 

Given,  now,  that  it  is  sufficient  to  work  on  scale  j  =  0  to  understand  exactly 
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what  the  wavelet  dmvative  does,  one  must  understand  the  ramifications  of  applying 
the  matrix  /IP  to  the  vector  sq.  In  the  next  subsection  the  similarity  between  the 
above  defined  matrix  If*  and  finite  difference  formulas  for  taking  the  derivative  will 
be  explored. 

4.4  Wavelet  Derivatives  and  Finite  Difference 

As  the  previous  subsection  illustrated,  the  essential  properties  of  the  wavelet  deriva¬ 
tive  are  contained  in  the  elements  of  the  matrix  R.  Recall  that  R  is  the  matrix  form 
of  the  mapping  from  sq  to  sq.  The  surprising  property  that  the  matrix  R  exhibits  is, 
however,  that  it  can  also  differentiate  evenly-spaced  samples  of  a  function.  That  is, 
R  acts  as  a  finite-difference  operator  when  applied  to  the  samples  of  a  function. 

This  subsection  is  in  three  parts: 

1.  The  similarity  between  wavelet  derivative  coefficients  and  finite  difference  coef¬ 
ficients  is  noted. 

2.  The  finite  difference  accuracy  of  the  coefficients  {rj}  derived  in  [9]  will  be  illus¬ 
trated,  and  it  will  be  proved  in  general  that  the  coefficients  {n}  can  differentiate 
polynomials  exactly  up  through  order  2Af  for  coefficients  {r/}  that  were  derived 
for  Daubechies  wavelets  D^m- 

3.  In  the  finite  element  method  under  certain  conditions  one  achieves  a  very  high 
order  of  accuracy  called  ’superconvergence.’  In  wavelet  differentiation  a  similar 
phenomenon  is  encountered.  This  phenomenon  is  defined  and  a  short  explana¬ 
tion  is  offered. 

4.4.1  Finite  Difference  Coefficients 

First  of  all,  it  is  useful  to  simply  note  the  similarity  between  the  coefficients  of  centered 
finite  difference  formulas  and  the  coefficients  used  to  construct  the  matrix  R.  The 
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Ord«r  oi  Accuracy 
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Table  2:  Optimal  centered  finite  difference  coefficients  with  order  of  accuracy. 
Wavelet  Convolution  Coefficients 
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Table  3:  Scaling  function  derivative  convolution  coefficients  for  Daubechies  wavelets. 


following  is  a  table  of  centered  finite  difference  coefficients  and  the  order  of  accuracy 
of  the  approximation  to  the  derivative: 

Recall  that  the  elements  of  the  matrix  R  derived  in  [9]  provide  the  transformation 
from  scaling  function  coefficients  of  a  function  to  the  scaling  function  coefficients  of 
the  derivative  of  the  same  function.  The  elements  of  R  for  the  and  D4  wavelet 
derivatives  are,  as  is  shown  in  the  following  table,  exactly  the  same  as  the  coefficients 
for  the  2-nd  and  4-th  order  centered  finite  difference  formulas.  Note  that  the  wavelet 
filters  become  quite  long  with  increasing  order.  Therefore,  only  the  right  side  of  the 
filter  will  be  shown  keeping  in  mind  that  these  filters  are  antisymmetric: 

The  fractions  for  the  D%  and  D%  wavelets  are  exact  but  complicated  and  provide 
little  insight.  Compare  the  following  decimal  representations  of  the  6-th  and  8-th 
order  finite  difference  operators  to  the  decimal  representations  of  De  and  Dg.  Once 
again,  only  the  right-hand  side  of  these  antisymmetric  filters  is  shown: 


The  coefficients  for  the  Dg  and  D9  derivatives  are  not  the  same  as  the  coefficients 
for  the  6-th  order  and  8-th  order  centered  finite  difference  derivatives,  but  the  dif¬ 
ferences  are  not  large.  Surprisingly,  however,  D9  has  the  same  accuracy  as  the  6-th 
order  finite  difference  operator,  and  Og  has  the  same  accuracy  as  the  8-th  order  finite 
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FD-6  and  D«  CoefficieoU 

FD-6  0  .750  -.150  .017 

Dt  0  .745  -.145  .015  .0003 

T»ble  4:  Comparison  between  numerical  values  of  optimal  6th  order  centered  finite 
difference  co^cients  and  Daubechies  6  scaling  function  derivative  convolution  coef¬ 
ficients. 


FI>>8  and  D%  Coefficients 

FD-8  0  .80  -.20  .038  -.0036 

Da  0  .79  -.19  .034  -.0022  -.0002  .0000008 

Table  5:  Comparison  between  numerical  values  of  optimal  8th  order  centered  finite 
difference  coefficients  and  Daubechies  8  scaling  function  derivative  convolution  coef¬ 
ficients. 

difference  operator.  §4.4.2  will  establish  this  accuracy. 

4.4.2  Finite  Difference  Accuracy 

To  establish  the  finite-difference  accuracy  of  the  wavelet-based  differentiation  coeffi¬ 
cients  note  that  a  centered-finite-difference  derivative  approximation  with  2K  anti¬ 
symmetric,  r*  =  — r_fc  implying  ro  =  0,  coefficients,  (rk)ks-Ky  ^  written 

-/>-.)•  (86) 

If  the  above  equation  is  exact  for  f{x)  =  for  9  =  0, ...,  N  but  not  for  9  =  ^  -f-  1 
then  the  equation  is  said  to  be  N-th  order  accurate.  Therefore,  one  must  check  to  see 
if 

9*r^  =  1C  »•*(*’+*  -  (87) 

when  f(x)  =  x^.  To  simplify,  one  can  let  Xj  =  j  and  check  the  following: 

=  E  ’■>(0  +  It)’  -  U  -  *)')•  (88) 

fc=l 

Now,  treating  the  coefficients  derived  in  [9]  as  nothing  more  than  finite-difference 
coefficients  one  can  check  the  accuracy.  The  following  table  contains  the  results  of 
applying  the  coefficients  from  [9]  to  various  polynomials; 
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Tkble  6:  The  degree  of  polyiMMnieb  differentieted  exactly  by  varknu  Daubechiee 
scaling  fiinctkm  derivative  coefficients. 

The  pattern  in  this  table  is  obvious  and  leads  to  the  fd^owing  theorem: 
Theorem:  If  ^z)  is  the  scaling  function  for  the  Daubechies  wavelet  denoted 
Iqr  ,  where  M  is  the  number  dt  vanishing  moments  of  the  wavdet,  then  the 
coefficients  {rj  derived  from 


and  Implied  the  to  evenly-spaced  samples  of  a  function  act  as  a  finite  difference 
derivative  operatm’  of  order  2M. 


The  proof  of  this  theorem  requires  two  resrilts,  as  well  as  the  Fourier  TVansform 
of  ^x).  First  will  be  found  followed  by  the  two  results  which  are  needed  and 
stated  in  themem  form.  Perhaps  the  first  proof  would  suffice  assuming  the  second 
result  is  weU-known.  However,  to  be  complete  the  second  result  is  also  proved. 

Ehurier  IVaiurfbrm  of  ^(x):  Recall, 

^x)  =  V2  52  hk^2x  -  k). 

Define  the  Fourier  Transform  of  ^(x)  as, 

ko = r 

^—00 
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TImmCopB) 

lt~l  MO 

^2*  -  *)«“*<*»■ 

teO  •'"* 

Lirt  p>m2»-’k  which  impliec  dx  s  (fy/2  to  get 

or  nmply 

=  d(f/2)4«/2). 

recalHngthat  ^(0  =  ^Efco  FVirthcnnore,  weget  =  W((/2)J7(^/4)^((/4). 

This  implies, 

but  ^*)  is  iKHrmalized,  ^0)  =  f  ^z)dx  ss  1.  Therefore, 


Theorem:  The  Fourier  Transform  of  {rj}  is  of  the  form, 

i^iere  c  6  C7  is  some  constant,  and  ’h.o.t.’  denotes  higher-order  terms  and  will  be 
used  again  in  the  proof. 

Proof:  B^n  with  the  expression  for  {n}: 

1 

If  we  <tefine, 

/(*)  = 
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and 


/Ky)  =  /*^(y) 

where  is  the  coovolution  <^>erator.  The  convolution  theorem  states  that  the  Fourier 
IVansform  (continuous  or  discrete)  of  a  convolution  is  the  product  of  the  Fourier 


TVansforms: 


If  we  define  ri  as, 


m = himy 


ri  =  p(/) 


then  the  semi-discrete  Fourier  transform  of  n  is 


where, 


r(0=  £  PU  +  2irfc). 

Ins— oo 


P(0  =  r  pix)e^dx, 

OO 


=  E 

ik=— oo 

Calculate  the  needed  Fourier  IVansforms  to  get. 


m = ko. 


where  ~  denotes  conjugation,  and 


m = wc). 


Combine  these  results  to  get  the  Fourier  IVansform  of  {pi}: 


m  =  \koH^ 
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Now»  we  need  to  know  tlM  brii»vi<Hir  <rf  Recall  frwn  the  definitimis  that, 

HU)  =  (1  + 

whexe  Q(e^)  does  not  have  poles  or  zeros  at  ^  =  x,  see  [2J.  That  is,  H{()  has  a  zero 
of  wder  Af  at  ^  =  ir.  Therefore,  H{(  +  *■)  has  a  zcto  of  order  Af  at  ^  =  0,  i.e., 

HU  +  ir)-c(^  +  h.o.t., 

then 

+  ir)p  =  +  h.o.t., 

where  a  s:  |c|^.  Recall  from  the  definitions  that, 

|f^(0P  +  (^«  +  ^)P  =  i- 

Combine  the  two  previous  relations  to  get, 

That  is, 

^|H({)l»U  =  0, 

for  n  =  1, 2Af  —  1.  The  Fourier  IVansfonn  of  <f>{x)  was  found  above: 

k()=n^(§)- 

We  get  an  expression  for  |^(0P  fro®, 

or 

itoi’  =  ft  iH(h’- 

i=i  ^ 

Now,  derivatives  of  this  expression  have  the  form, 
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and  tme  can  see  that  if,  ^  —  1,  then 

=  0, 

fcNT  n  =  l,...,2Af  —  1.  FVom  this  information  we  can  see  that  a  series  expansion  of 
1^01^  about  (  =  0  would  be  of  the  form, 

WOI’  =  «  +  »("'  +  k.o.t. 


But,  ^(x)  is  normalized  implying  that  ^0)  =  1  and  therefore  |^(0)p  =  1.  The 


expansion  becomes. 


itoi’  =  1 + ‘f’" + 


where  6  6  C.  Recall  that  we  are  looking  for  the  semi-discrete  Fourier  transform  of 
f(^),  which  we  see  from  above  is, 

fK)=  E  (>«  +  2>rt)=  f:  i(f  +  +  2it)|». 

*=—00  *=— oo 

We  now  need  to  find  the  behaviour  of  +  2xfc)p  when  k  ^  0.  Note  that  in  the 


expression. 


the  arguement. 


i^(f+2xt)i'=nifl(i^)i>. 

i=i  ^ 


will  for  some  j  be  equal  to 


^  -l-2xfc  _  j1  , 

^  2^-^  ’ 


(  +  2Kk_  ( 


modulus  2x.  That  is,  if  k  is  odd  in  the  expression  (kx)/(2^~^)  then  stop  when  j  =  1 

A  A 

since  we  can  subtract  multiples  of  2x  from  (kir)  without  changing  H  since  H  is  2x 
periodic.  If  k  is  even,  then  for  some  j,  the  number  (fc)/(2^“^)  will  be  odd  at  which 
point  we  ag»n  subtract  some  multiple  of  27r.  Consequently,  in  the  infinite  product, 

J=1  ^ 
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ihae  will  always  be  a  term,  when  ib  0,  on  the  right  hand  side  with  the  f(Nrm, 
\H{^  +  *•)!’.  But,  from  above,  we  found  that, 

|H{|  +  »)l’  =  OK”*)  + 

But  this  implies  that  for  Jb  ^  0  that  the  contributions  to  the  infinite  sum, 

f;  !*({+ 2x1)1’ 

ks-oo 

are  of  That  is, 

f;  \^i(  +  2irk)\^  =  l-^b(^^  +  h.o.t. 

k=—oo 

Ultimately,  we  need  the  semi-discrete  Fourier  transform  of  {r}: 

^(0  =  *  S  (^  +  2irfc)|<^(^  -I-  21rl;)|^ 

ib=-oo 

or 

H()  =  ‘(  E  |^((  +  2xt)|’  +  2xi  f;  t|^(f  +  2x*)|’. 

*=—00  *s— 00 

We  already  know  the  behaviour  of  the  first  term  on  the  right-hand  side.  The  second 
term  on  the  right-hand  side  cannot  contribute  powers  of  ^  lower  than  2M  since  it 
differs  from  the  summation  in  the  first  term  only  by  a  multiple  of  k  which  does  not 
allow  the  low  power  contribution  when  k  =  0.  The  final  step  is  to  illustate  the  second 
term  on  the  right-hand  side  is  an  odd  function  implying  that  the  lowest  power  of  ^  it 
can  contribute  is  2M  -f  1,  the  first  odd  number  past  2M.  That  is,  define 

/«,t)=t|^(<  +  2xt)|’, 

and  note  that  in  the  infinite  summation  that  there  is  always  a  term  with  positive 
k  and  a  term  with  negative  k.  The  summation  of  all  such  -f-A;  terms  and  —k  terms 
yields  odd  functions. 

m,  k) + /u,  -ib) = /{-e,  k)  ±  /(-^,  -k), 
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and  this  implies  the  desired  result  leaving  us  with  the  conclusion  of  the  proof  that, 

f{0  =  ii  +  ibe^*'  +  h.o.t. 

This  completes  the  proof.// 


Lemma:  Let  {rj}  be  a  finite  set  of  coefficients.  These  coefficients  can  be  called 
the  coefficients  of  a  finite  difference  approximation  to  a  first  derivative,  or  these 
co^cients  can  be  called  a  finite  impulse  response  filter,  or  FIR  filter.  FWthermore, 
let  the  coefficients  be  antisymmetric:  n  =  — r_{  which  implies  ro  =  0.  If  the  Discrete 
Fourier  Transform,  or  DFT,  of  {n}  is  of  the  form 

r(0  =  *^  +  cr'^*  +  ^-‘>.t.,  (89) 

for  some  constant  c  €  C,  then  the  filter  {r/}  when  applied  to  evenly-spaced  samples 
of  a  function  can  differentiate  in  a  finite  difference  sense  with  accuracy  of  order  m. 
That  is,  {r/}  can  differentiate  polynomials  exactly  up  to  a:”*. 

Before  the  proof,  note  that  the  DFT  of  a  filter  which  is  designed  to  approximate 
differentiation  in  the  space  domain  should  approximate  in  the  frequency  domain: 

(90) 

ax 

That  is,  differentiation  filters  are  attempting  to  approximate  the  frequency  of  a  pure 
sinusoidal  mode. 

Proof: 

Let  the  DFT  for  {r/}  be  defined  as, 

HO  =  E’-i'*"-  (91) 

I 

Break  up  the  summation  to  write  as, 

HO  =  ^0  +  H  (92) 

J>1  /<-! 


39 


Now,  impoae  the  antisymmetry  to  get, 


l>i 


(93) 


Using  the  series  expansion  about  zero  for  the  complex  exponential  one  gets, 


or  factor  to  get. 


1>1  k=o 

Interchange  the  summations  to  get, 


i>i  k=o 


r(0  =  E'-,f;l^(/* -(-/)*). 

IX.S 


(94) 


(95) 


‘‘(«=f:^E’-‘('‘ -(-')*)• 

kso  *•  I>1 


(96) 


Recall  that  the  hypothesis  was  that  r(^)  =  if  +  cf™'*'^  +  h.o.t.  This  implies  that. 


o  =  En(<‘ -(-')*)  (97) 

<>» 

for  A;  =  0  and  for  2  <  k  <m.  Furthermore,  1  =  J^/>i  r/(i*  —  (—/)*)  must  hold  when 
k  =  1.  But  these  are  exactly  the  conditions  that  must  hold  for  a  filter  to  be  able 
to  differentiate  ex2u:tly  polynomials  up  through  order  m  which  are  centered  at  zero. 
The  proof  for  polynomials  centered  at  some  arbitrary  position  requires  a  shift  in  the 
index  but  the  results  are  the  same.  This  completes  the  proof./ / 


This  completes  the  theorems  which  are  at  the  heart  of  the  paper.  The  next  section 
discusses  the  high  order  of  accuracy  of  the  coefficients  {rj. 

4.4.3  *Superconvergence’ 

Note  that  the  matrix  can  differentiate  exactly  polynomials  up  to  degree  2M  for  the 
Daubechies  wavelet  D2M  when  thought  of  as  a  finite-difference  operator,  even  though 
the  scaling  function  subspace  Vo  can  only  represent  exactly  polynomials  only  up  to 
degree  Af  —  1.  A  Similar  phenomenon  is  encoimtered  in  the  finite  element  method 
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under  partictilar  choices  of  the  approximation  grid  and  is  known  as  superconverg^ce 

111]. 

To  understamd  the  source  of  superconvergence  in  the  wavelet  derivative  it  is  suf¬ 
ficient  to  have  a  good  imderstanding  of  the  proofs  in  the  previous  subsection.  Let  us 
note  the  sources  of  the  powers  of  (  in  the  expression  for  the  DFT  of  the  coefficients 

{n}: 

^(0  =  *^  +  -f  h.o.t. 

Recall  the  definition  of  {rj}, 

n  =  /  <f>ix  -  l)^(f>{x)dx, 

•/— oo  uX 

as  well  as  the  definitions  of  and  ^<f>{x)‘. 

“  0. 

{=0 

kx)  =  2  X)  ^(2®  -  0- 

1=0 

The  sources  of  the  powers  of  ^  are  now  apparent:  M  powers  come  from  and  M+ 1 
powers  come  from  The  ’superconvergence’  for  the  wavelet  derivative  can  be 

explained  by  the  similarity  between  the  equations  which  define  ^(x)  and  i<f>{x).  That 
is,  they  are  defined  by  dilation  equations  which  differ  only  by  a  multiple  of  2. 

The  next  section  will  smnmarize  and  conclude  this  paper. 
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5  Conclusion 


A  restatenKait  of  the  thesis  of  this  paper  is  given  first  followed  by  a  brief  outline  of 
the  argument. 

Given  the  evenly-spaced  samples  of  a  {>eriodic  function,  /,  then  the  matrix 
derived  for  the  Daubechies  wavelet  has  the  effect,  when  applied  to  /,  of  a  finite- 
difference  derivative  operator  of  degree  2M. 

The  hear  of  the  argument  of  this  paper  is  contsuned  in  §(3)  and  §(4).  In  §(3)  it 
was  established  that  if  given  the  evenly-spaced  samples  of  a  periodic  function  /(x) 
then  the  scaling  function  coefficients  so  of  the  function  at  the  finest  scale  can  be 
approximated  by  a  quadrature  formula  which  in  matrix  form, 

do  =  Cf, 

yields  a  circulant  matrix  C,  where  do  approximates  sq.  F^irthermore,  in  §(3)  it  was 
noted  that  all  circulant  matrices  with  the  same  dimensions  commute.  In  §(4)  it  was 
noted  that  the  coefficients  which  map  the  scaling  function  coefficients  at  the  finest 
scale  of  a  periodic  fimction  to  the  scaling  function  coefficients  at  the  finest  scale  of  the 
derivative  of  the  function  is  also  circulant  in  form  when  written  in  matrix  notation, 

Po  = 

Furthermore,  it  was  observed  that  the  matrix  can  differentiate  evenly-spaced  sam¬ 
ples  of  a  polynomial  in  a  finite-difference  sense  exactly  up  to  the  order  of  the  wavelet. 
Also,  when  i?  is  applied  to  the  evenly-spaced  samples  of  a  periodic  function  then  If" 
is  circulant.  Now,  combine  the  results  of  §(3)  and  §(4)  to  get  the  following  relation: 

P  =  C-^FPCf. 

Throughout  the  paper  it  has  been  noted  that  C  and  BP  are  circulant  in  form  when 
f{x)  is  periodic  and  that  circulant  matrices  commute.  Therefore,  the  previous  relation 
simply  becomes. 
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where  BP  is  acting  as  a  finite  difference  operator. 

A  note  concerning  notation  is  in  order.  In  the  introduction  the  matrices  C,  D, 
and  the  differentiation  matrix  T>  were  defined.  Under  the  scenario  developed  in  this 
paper,  the  wavelet  matrix  C  is  the  same  as  the  matrix  C  from  the  introduction.  The 
matrix  D  from  the  introduction  becomes  the  matrix  BP.  Likewise,  the  matrix  V  is 
also  BP  since  for  wavelets  C  and  BP  commute.  That  is,  for  evenly-spaced  samples 
and  periodicity  W  is  the  wavelet  differentiation  matrix  which  has  the  effect  of  a  finite 
difference  operator. 

The  importance  of  the  thesis  of  this  paper  is  that  under  periodicity  and  an  evenly- 
spaced  grid  one  can  understand  the  wavelet  differentiation  matrix  in  terms  of  a  finite 
difference  operator  with  the  accuracy  given  by  the  superconvergence  theorem. 
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Appendix  A:  Wavelets  Supported  on  (0,3M) 

In  this  appendix  our  wavelets  are  supported  on  [0, where  M  is  the  number  of 
vanishing  moments  of  the  wavelet.  These  are  not  the  usual  Daubechies  wavelets,  but 
for  these  wavelets  the  scaling  function  coefficients  of  a  periodic  function  /(z)  can  be 
approximated  with  error  of  order  M  simply  by  sampling  /(z)  at  the  correct  location. 

To  begin,  assume  that  there  exist  a  unique  tm,  fixed  for  a  fixed  number  of  vanishing 
moments,  Af,  of  the  wavelet,  such  that 

j  4{x  +  TM)x”^dx  =  0 

for  m  =  1,2,  ...,Af  —  1.  Furthermore,  recall  the  definition  of  the  scaling  function 
coefficient  and  expand  the  integrand  /(z)  in  a  Taylor  series  about  zq  (/q  =  /'(zq)): 

=  J  /(x)^(z  -  k)dx  = 

fo  J  <(t(x  -  k)dx  +  /q  j{x-  Xo)<f>ix  -  k)dx  -f  fH  xo)V(®  -  +  .... 

Now,  shift  the  variable  of  integration  by  y  =  z  —  t  —  Ar,  and  choose  the  point  of 
expansion,  zq,  to  be  r  4- 1?  to  get, 

= 

/(r  +  fc)  y  <t^{y  +  T)dy  +  f\y  +  k)j  y^(y  +  T)dy  +  f{y  +  k)J  y  V(y  +  T)dy  +  .... 
Now,  rename  t  as  tm  and  the  above  integrals  are  of  the  form, 

J  <f>{x  +  TM)x”"dx  =  0, 

and  therefore  vanish  for  m  =  1, ...,  M  —  \  leading  to, 

=  f{TM  +  k)-\-  +  k)j  y"^(y  +  TM)dy  +  ..., 

i.e.,  the  approximation  of  the  scaling  function  coefficient  up  to  order  M  is  made 
by  sampling  /(z)  at  the  position  tm  +  k. 
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Note  thet  ell  of  the  above  cakuUtkMU  could  have  been  carried  out  for  the  first 
dorivative  of  /(«)  giving  an  approrimation  to  the  scaling  function  co^cients,  of 

r(*): 

i*  =  /(r  +  k)  +  +  fc)  y  y"^(y  +  T)dy  +  .... 

It  was  assumed  above  that  there  exist  one  rv  such  that 

y  +  m)x”"dx  =  0, 

for  m  =  1, ...,  A#  —  1.  For  m  =  1  this  tji/  is  easy  to  find: 

j  4(x  +  TM)xdx  =  j  ^(x)(x  -  rsi)dx 

=  j  x<ft{x)dx  J  ^(x)rfx. 

But  f  <^{x)dx  =  1,  therefore, 

t-m  =  y  x^(x)dx. 

That  is,  Tftf  is  simply  the  first  moment  of  ^(x).  To  find  tm  for  m  >  1  the  calculations 
are  simple  but  a  bit  longer  and  require  the  result  from  the  following  theorem  to  show 
that  there  is  one  tm  which  is  the  same  for  all  m  =  1, ...,  Af  —  1. 

ltJi>ix)dx  =  1  and  there  exists  T  such  that /^(x+T)x"*dx  =  Oform  =  l,...,Af— 1 
then  /^(x)x”*dx  =  (/^(x)xdx)”*  for  m  =  1,...,  Af  —  1. 

Proof:  Start  with 

y  <ff{x  +  T)x”*dx  =  0, 

and  let  y  =  X  +  r  to  get, 

/  «>(!/)(»  -  t)”  =  0. 

Using  the  binomial  theorem  this  becomes, 

/«»)£(  7 
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Let  the  moments  of  ^(z)  be  denoted  by  Mi  =  f  ^(x)z'dz  to  get 

E  ( 7 )  (-’■)”"*''  =  “• 

A  simple  calctilation  yields  t  =  Mi.  Using  this  value  of  t  and  summing  only  up  to 
m  —  1  the  previous  expression  becomes, 

m— 1 


E  ( 7 )  = 0. 


Or, 


«» =  -  E  ( 7 )  (-ir-'wr-'M,. 

Prom  the  hypotheses  it  is  known  that  Mo  =  f  ^(z)dz  =  1.  Therefore,  Mp  =  Mi  for 
p  =  0, 1,  and  with  this  knowledge  it  is  easy  to  show  that  Mp  =  Mi  for  p  =  2: 

=  -  E  ( 7 ) 

which  holds  for  m  =  1,2.  Combine  the  powers  of  Mi  to  get, 

=  -«r  E  ( 7 ) 

But,  this  is  nothing  more  than. 

Mm  =  -Mn(i  -  ir  - 1], 


or  simply. 


Mm  =  Mr, 

where  m  =  1,2.  The  proof  is  complete,  since  higher  powers  of  m  can  be  found  by 
induction. 
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Af^peodix  B:  Moments  of  the  Scaling  FHmction 

In  this  appendix  the  m<»nent8  of  ^(x)  will  be  calculated  in  closed  fonn.  B^n 
with  the  definition  of  the  scaling  function, 


^x)  =  5^^(2x  -  i:). 

k 

Next,  calculating  the  m-th  moment  of  ^(x)  yields. 


J  ^(x)x"*  =  53  ^*/  ^(2®  -  k)x”^dx. 


Change  the  variable  of  integratimi  such  that  y  —  2x  —  k  to  get, 


/  = E  **/  «»)(i/2r(» + 

=  (1/2)”*'  E*‘/  Miy  + 


Now,  recall  the  binomial  theorem  to  get. 


j  = (i/2r«  £4,  /  mt  ( 7 )  v'*”"'*'!' 

Rewrite  the  moments  of  ^(x)  as  M{  =  /  x^^(z)dx  to  get, 

Af«  =  (1/2)"*+*  E  f  7 ) 

1=0  \  /  fc 

Now  let  m  =  hkk^  to  get 

= (i/2r+'  E  ( 7 ) 

Now,  Mm  can  be  defined  in  terms  of  Mi  for  »  =  0, ...,  m  —  1: 


M„(2”«  -  2)  =  £  (  7  ) 

Note  that  the  moments /t/  can  be  found  by  direct  calculation  given  that  the  Daubechies 
filter  coefficients,  h/k,  are  already  known. 

The  moments  of  ^(x)  can  now  be  used  to  find  the  mapping  from  the  evenly- 
spaced  samples  of  a  function  /(x)  to  the  scaling  function  coefficients.  In  section  3 
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i  Mj  ftj 

0  1  2  .3661  . 

1  .63395  1.2679  .6340 


Table  7:  Scaling  function  moments  for  the  Daubechies  4  wavelet. 


i 

Mi 

t^i 

Ci 

0 

1 

2 

-.0235 

1 

1.0054 

2.0108 

1.0426  . 

2 

1.0108 

2.0216 

-.0199 

3 

.90736 

.5078 

.0009 

Table  8:  Scaling  function  moments  for  the  Daubechies  8  wavelet. 

this  mapping  was  denoted  in  matrix  form  as  the  matrix  C.  The  elements  c,  which 
define  the  rows  of  this  matrix  have  already  been  given  for  the  wavelet  Dq.  The 
comparable  results  for  the  D4  and  Dg  wavelets  are  given  in  the  accompanying  tables. 
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